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I’relace 


Under  Conlracl  No.  F49620-C-89-()038,  N'I  NF/NORSAR  i.s  conducting  research  within  a  wide 
range  of  subjects  relevant  to  seismic  monitoring.  The  emphasis  of  the  research  program  is  on 
developing  and  assessing  methods  for  proce.ssing  of  data  recorded  by  networks  of  small-aperture 
arrays  and  3-component  stations,  for  events  both  at  regional  and  tele.seismic  distances.  In  addition, 
more  general  seismological  research  topics  are  addressed. 

Each  quarterly  technical  report  under  this  contract  presents  one  or  .several  separate  investigations 
addressing  specific  problems  within  the  .scope  of  the  .statement  of  work.  Summaries  of  the  research 
efforts  within  the  program  as  a  whole  are  given  in  annual  technical  reports. 

This  Scientific  Report  No.  10  prc.senLs  a  manu.scripl  entitled  "Diffraction  and  seismic  tomography", 
by  D.J.  Doombos. 
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Summary 


Diffraction  tomography  is  formulated  in  such  a  way  that  the  data  (travel  time  -  or 
waveform  perturbations)  are  related  to  the  medium  perturbations  through  the  sum  of 
two  terms.  The  first  term  is  the  ray  integral  of  ordinary  tomography  and  involves  only 
phase  perturbations.  The  additional  diffraction  term  involves  both  phase-  and 
amplitude  perturbations.  The  diffraction  term  is  linear  in  the  gradients  of  the  velocity 
perturbation  in  an  acoustic  medium,  the  gradients  of  the  elastic  and  density 
perturbations  in  an  elastic  medium,  and  the  gradients  of  the  boundary  perturbations 
the  wave  is  crossing.  This  formulation  has  the  additional  advantage  that  unwanted 
diffractions  from  the  nonphysical  boundary  of  the  region  under  study  can  be  easily 
removed.  Acoustic  scattering,  elastic  scattering,  and  scattering  by  boundary 
perturbations  are  analysed  separately.  Attention  is  paid  to  the  adequacy  of  the 
acoustic  approximation,  and  to  the  difference  between  perturbations  of  a  boundary 
level  (topography)  and  perturbations  of  boundary  conditions.  These  differences  are 
irrelevant  for  ordinary  seismic  tomography.  All  results  are  based  on  first-order 
approximations  (Born  or  Rytov),  as  is  the  case  for  other  published  methods  of 
diffraction  tomography. 
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Introduction 


Tomographic  methods  have  played  an  increasingly  important  role  in  many  fields 
where  ray  theory  is  applicable  to  describe  the  propagation  of  waves.  Very  efficient 
inversion  algorithms  have  been  developed  for  situations  where  sources  and  receivers 
are  placed  at  regular  intervals  in  a  homogeneous  background  medium,  as  is  common 
in  X-ray  medical  applications  (Barrett,  1984).  Similar  algorithms  have  been  used  in 
seismic  cross-well  tomography  (see  e.g.  Worthington,  1984,  and  references  therein). 
However  in  seismological  applications,  sources  and  receivers  are  usually  located  at 
irregular  intervals,  and  the  homogeneous  background  medium  is  often  not  a  good 
assumption.  To  some  extent  the  same  may  be  true  for  cross-well  tomography.  Thus 
seismic  tomography  usually  cannot  take  direct  advantage  of  fast  transform  methods. 
However  the  problem  of  travel  time  inversion  can  be  posed  in  the  form  of  a  general 
linearized  inverse  problem,  and  in  this  form  seismic  tomography  has  been  widely 
applied  both  for  exploration  purposes  (Worthington,  1984;  Ivansson,  1987),  in 
lithospheric  studies  (Aki,  Christoffersson  and  Husebye,  1977)  and  to  study  the  Earth’s 
deep  interior  (Dziewonski,  Hager  and  O’Connell,  1977). 

Despite  the  widespread  application  of  these  methods,  a  number  of  limitations  were 
recognized  at  an  early  stage.  A  fundamental  limitation  is  due  to  the  neglect  of  wave 
diffraction  off  the  geometrical  ray  path.  This  neglect  limits  the  resolution  of 
tomographic  results  in  general.  The  limitation  is  especially  serious  in  circumstances 
where  diffraction  effects  dominate  the  observations,  as  may  be  expected  in  the 
presence  of  low-velocity  regions  among  others  (Wielandt,  1987).  To  overcome  these 
limitations,  diffraction  effects  have  been  taken  into  account,  in  a  first-order 
approximation,  in  methods  called  diffraction  tomography  (e.g.  Devaney,  1984).  It  was 
later  recognized  that  these  methods,  when  applied  to  reflection  data,  are  closely 
connected  to  wavefront  migration  (Miller,  Oristaglio  and  Beylkin,  1987).  However  in 
accord  with  the  original  formulation  of  tomography,  current  methods  in  diffraction 
tomography  assume  a  homogeneous  background  medium  and  regular  sampling  of  the 
wave  field,  such  that  Fourier  transform  methods  can  be  used  for  the  solution  of  the 
scattering  problem  at  hand.  In  its  present  form  these  methods  are  of  limited  use,  at 
least  in  earthquake  seismology. 

In  this  paper,  diffraction  tomography  will  be  formulated  in  such  a  way  that  the 
observation  variable  (a  travel  time  residual,  or  more  general  a  waveform 
perturbation)  is  approximated  by  the  sum  of  a  geometrical  ray  term  which  itself  is  the 
basis  of  ordinary  seismic  tomography,  and  a  diffraction  term  which  is  an  integral  over 
the  volume  of  heterogeneity  and  over  the  surface  of  boundary  perturbations.  The 
integrand  is  linear  in  the  gradients  of  heterogeneity  perturbations  (velocity,  elastic 
constants,  and  density)  and  the  gradients  of  the  boundary  perturbations  sampled  by 
the  wave  field.  This  means  that  diffraction  tomography  can  still  be  posed  in  the  form 
of  a  linearized  inverse  problem.  TTie  results  are  based  on  first-order  scattering  theory 
as  is  the  case  for  other  published  methods  of  diffraction  tomography.  In  the  following 
we  separately  consider  acoustic  and  elastic  scattering,  and  scattering  by  boundary 
perturbations. 
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Acoustic  scattering 


Consider  the  wave  motion  u(x,t)  in  a  medium  with  velocity  c,  slowness  s=l/c 
and  density  p.  The  reference  medium  has  velocity  c°,  slowness  s9=llc^,  and  the 
wave  motion  is  u°(x,t).  The  slowness  perturbation  in  a  volume  V  is  Ss=s-sf^,  and 
the  ensuing  scattered  wave  motion  is  Su=u-u^.  No  density  perturbations  are  assumed 
here;  acoustic  scattering  with  density  variations  has  been  treated  by  Stolt  and  Weglein 
(1985);  perturbations  in  both  the  density  and  in  the  elastic  parameters  are  treated  in 
the  next  section,  on  elastic  scattering.  A  standard  procedure  using  the  Born 
approximation  gives  (eg.  Aki  and  Richards,  1980): 


6 


u(x,  t)  -- 


I 


2pc  fis  G*U°  dv 


(1) 


where  G  is  the  Green’s  tensor  for  the  medium.  Let  u®  and  G  be  given  by  ray 
theory: 


t)  =  vVC)  f(t-T°)  (2) 

G{x,^,t)  ==v(x)vi'(5)  AHx,^)  b(t-T^)  (3) 


Quantities  as.sociated  with  the  incident  and  scattered  wave  have  superscripts  °  and  \ 
respectively;  v,  v®,  v'  are  unit  displacement  vectors,  and  /I®,  A'  are  amplitude 
factors.  A°  and  A^  are  allowed  to  be  complex,  but  this  requires  a  generalization  of 
equations  (2)  and  (3)  such  that  for  example  A^(t)  is  replaced  by  Re(A^)f(t)  - 
Im(A^)g(t),  where  g(t)  is  the  Hilbert  transform  of  f(t).  We  assume  this 
generalization  in  all  the  following  equations  if  required,  although  for  simplicity  we 
have  not  brought  this  out  in  the  notation.  Substitution  in  equation  (1)  gives 


5u(x,  t)  - 


^2pc  (x) 


(v*.v»)  5s  dv 


(4) 


which  is  a  well-established  result  (c.f.  Coates  and  Chapman,  1990). 
F^'ach  component  of  Su  can  be  written  in  the  form 


where 


r  =  P+r 


and  V  =  2pc  A°A^  v^(x) 


(6) 


Equation  (5)  is  transformed  into 


6u^  (x,  t) 


- 1  £(t-x)  J  w  (V®.V^) 


6s  J  dS^dx 


(7) 


where  S^.  is  the  surface  T=constant,  and  J  is  the  Jacobian  of  the  transformation 
(see  Fig.  1).  The  transformation  to  an  integration  over  what  is  called  isochron 
surfaces  5^  has  been  applied  before  (e.g.  Haddon  and  Buchen,  1981;  Miller, 
Oristaglio  and  Bcylkin,  1987;  Cao  and  Kennett,  1989).  The  lower  integration  limit  in 
the  T  integral  is  the  stationary  travel  time  of  the  geometric  ray,  and  we  have  chosen 
the  upper  limit  such  that  the  bounding  surface  of  V  coincides  with  for  constant 
r„.  In  the  notation  it  is  implicit  that  we  have  assumed  the  stationary  travel  time 
to  be  a  minimum  with  respect  to  path  variations  due  to  scattering  points  in  V.  If 
is  a  maximum  with  respect  to  these  variations  we  have  to  reverse  the  integration 
limits,  and  for  a  "minimax"  time  (i.e.  is  a  saddle  point  with  respect  to  the  relevant 
path  variations)  the  integral  is  to  be  split,  but  the  final  results  will  be  essentially  the 
same.  Integration  by  parts  gives 

6Ui  {X,  t)  “  -  f  (t-Tj  /  W  (v«.v*)  6s  d 
+  /  rv  (v®.v^)  6sj^ds^.^ 

^xu 

-f  V  (v®.vM  f(t-T)  dv  (8) 


where  we  have  neglected  the  variation  of  W,  since  it  is  coupled  to  Ss  and  would 
produce  a  second  order  effect.  (However,  variations  in  W  would  have  to  be  retained 
in  a  similar  integral  equation  for  u®).  The  integration  surface  in  equation  (8) 
encloses  the  infinitesimal  volume  SV^  about  the  geometric  ray.  Let  the  ray  path  be 
parameterized  by  a  ,  with  limits  and  a,.  TTien 


da 

Oo 


where  6S^  is  an  infinitesimal  element  of  surface  normal  to  the  ray  in  a  (see  Fig.  1). 
Within  this  surface  we  can  expand  the  time  (Farra  and  Madariaga,  1987): 


^  ^  ^  ^  fig'  Cfiff 


(9) 


where  Sq  denotes  the  position  (in  ray  centred  c(X)rdinates)  within  65^,  and 


C  =  ‘  ‘ 


The  2x2  matrices 


3g(Q) 

5p(Oo) 


and  0^ 


(  dqLal 

\  dpiOj^) 


define  the  geometrical  spreading  between  Oq  and  a  and  between  a,  and  a, 
respectively  (Cerveny,  1985).  Here  p(a)  is  the  2-D  slowness  vector,  in  ray  centred 
coordinates.  For  minimum  travel  time  t„.  the  quadratic  form  in  equation  (9)  is 
positive  for  all  Sq,  and 


ftg*"  C  5g  =  a  =  constant 


is  the  equation  of  an  ellipse  with  area  7ro/|C|*/^.  Hence  the  surface  area 
bounded  by  t  =constant  can  be  obtained  by  putting  a=2(T-T ^=2Sr: 


6  S^(a)  =  2nbx/\C\^^^ 


(10) 


Using  this  result  in  equation  (8),  we  have 

■S’tn  =  -^1/2 


4 


Moreover,  following  Coates  and  Chapman  (1990): 


(12) 


where 


dq(a,)  ' 

^  ~  i5p(Oo)  ^ 

is  the  spreading  matrix  between  and  o^.  The  amplitude  factor  of  the  incident 
wave  (equation  2)  is 


A°  ~  ±  (pc|(?®|) 


(13) 


where  the  ±  sign  corresponds  to  the  minimum/maximum  time  character  of  the  ray 
between  and  a;  for  a  minimax  ray,  |g°|  is  negative  and  will  be  imaginary. 
The  result  (13)  needs  to  be  generalized  if  multiple  caustics  exist.  These  can  be  taken 
into  account  by  introducing  the  so-called  KMAH  index  a  (Chapman,  1985;  Coates 
and  Chapman,  1990)  and  replacing  /I®  by  \A^\  exp  (iual2).  The  index  a  increases 
by  an  integer  (normally  1)  each  time  the  ray  touches  a  caustic  surface. 

The  amplitude  factor  of  the  Green’s  function  (equation  3)  is 

(pcp^c^  \0^\)  (14) 

411 

The  rules  for  the  sign,  and  the  required  generalization  if  multiple  caustics  exist,  are 
the  same  as  for  /1°. 

Combining  equations  (6)  and  (11)-(14)  we  can  rewrite  the  first  term  of  equation  (8): 


5 


6u“(x,  t)  =  Tf  (C-tJ  j  W  (v®.v')  6s  d5^„ 

=  -  u“(x,  t)  j  bs  da 

Oo 

(15) 

where 

u^(x.c)  =Vj(x)  iA(x) 

(16) 

and  A(x)  -±  (PjCj  |0|) 

To  proceed  with  the  other  terms  of  equation  (8),  we  attach  to  each  point  in  K  a 

unit  vector  rj  normal  to  the  surface  j=constant  in  the  direction  of  increasing  t. 

For  ^  not  close  to  the  geometric  ray: 

T1  =  {p®-p^)  /  |P®-P^| 

(17) 

and 

dx/dr\  Ip®  pM 

(18) 

where  and  /»'  are  the  3-D  slowness  vectors  of  the  incident  and  scattered 

respectively.  In  the  second  term  of  equation  (8)  we  then  have 

waves, 

=  (dx^ldn)  =  |p*-p‘|'^ 

(19) 

In  the  third  tern  of  equation  (8): 

T1.V6S  /  =  (p“-pM  .V6s  /  1p®-pM" 

ox  dr) 

(20) 

Although  this  result  is  not  valid  close  to  the  ray,  it  is  obvious  from  symmetry 
considerations  that,  provided  S  s  is  continuous  in  any  small  surface  area  65, 
normal  to  the  ray  in  a  and  hounded  by  ST—comtanl: 

fa 

n 

dS  -  0 


/ 

at 


hence  vSs  does  not  contribute  in  a  region  bounded  by  Sr-conslant  close  to  the 
ray,  and  this  region  can  be  excluded  from  integration. 

Summarizing,  using  equation  (15),  (19)  and  (20)  we  rewrite  equation  (8): 


5Uj  (x,  t)-  -  u  "(x,  t)  j  bs 


da 


+  f  (  t-T„)  / 


tv  (v^.yt) 


_ ds 


-  f  h,  (V«.v‘)  f  (t-t)  dv 


(21) 


'Fhe  first  term  is  just  the  first  term  of  a  Taylor  series  expansion  of  u^(x,t)  due  to  a 
change  in  geometrical  travel  time.  Coates  and  Chapman  (1990)  obtained  the 
equivalent  pha.se  delay  term  in  the  frequency  domain,  by  a  different  method.  The 
second  term  expresses  diffraction  from  the  boundary  of  K  If  the  boundary  is 
nonphysical,  this  term  should  be  deleted  (the  implicit  assumption  being  fis=constant 
outside  C).  The  third  term  expresses  diffraction  due  to  changes  in  the  velocity 
perturbation  vSs.  Retaining  the  first  and  the  third  term,  (5u  ™  and  Su’*,  we  have  in 
the  frequency  domain 


6u.{x,(t>)  “  6C/"(x,  u)  +  6{//(x,  0)) 

=  10)  {/"(x,  0))  j  bs  do 

+  io)  i-’lo))  /  (V(v®.v^)  ^  exp(io)T)dV  (22) 

V  Ip^-pN^ 


The  Rytov  approximation  is  (Tarantola,  1987,  page  484); 


7 


In 


t7^(x,  0>) 


t//’(x,w) 


(x)  A(x) 


t//’(x,  w) 

^  W(v°.v^) 


i(l) 


/ « 


s  do 


(P  exp{iG)(t-Tj}  dl^ 


lpO-pl|2 


(23) 


The  first  term  in  the  brackets  is  the  usual  ray  integral  of  seismic  tomography.  The 
second  term  is  the  basis  of  diffraction  tomography.  The  integration  volume  V  can 
be  chosen  so  as  to  include  (an  appropriate  fraction  of)  the  Fresnel  zone.  We 
anticipate  that  the  degree  of  improvement  of  the  tomographic  image  depends 
primarily  on  an  accurate  estimate  of  the  phase  6>t  in  the  diffraction  integral.  Thus 
T  may  have  to  be  calculated  iteratively  using  previous  tomographic  results. 
Preliminary  results  suggest  that  inclusion  of  the  diffraction  integral  can  be  a  significant 
improvement  even  when  inverting  only  travel  time  residuals  (c.f.  Whitten  and  King, 
1990).  In  the  latter  case  only  the  real  part  of  the  integral  is  used  in  the  inversion;  the 
imaginary  part  gives  the  amplitude  perturbation. 

Any  of  the  velocity-  or  slowness  parameteriziUion  schemes  in  common  use  will  convert 
equation  (23),  (22)  or  (21)  to  a  form  that  is  linear  in  the  slowness  parameters,  so  that 
standard  inverse  methods  can  be  applied.  The  common  case  of  a  block 
parameterization  is  special  in  the  sense  that  the  diffraction  term  contributes  only 
diffractions  from  the  block  boundaries.  Thus,  at  the  boundary  5-,^  between  two 
adjacent  blocks  j  and  k: 


6s  =  -  (6s]  ^ 


where  is  the  surface  normal  pointing  into  block  k,  and  65^,,  is  the  delta- 
function  on  The  diffraction  term  is  then 


huf{x,  t)  =  S  ^  . 

J  k>j 


[6  s]  ^ 


/  W  (v“.vM  f(t-x)  dS 


)K 


(24) 


where  the  sum  over  k  is  restricted  to  blocks  adjacent  to  j,  and  the  additional 
restriction  k>j  is  used  to  avoid  duplicating  interfaces. 
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Elastic  scattering 

The  problem  of  first-order  elastic  scattering  (Born  approximation)  has  been  treated  by 
many  authors  (e.g.  Hudson,  1977).  Our  purpose  here  is  to  emphasize  the  agreements 
and  differences  between  the  results  for  the  elastic  and  the  acoustic  case.  In  particular, 
we  are  interested  in  the  range  of  validity  of  the  acoustic  scattering  assumption.  We 
assume  an  isotropic  medium,  and  use  6p,  Sk  and  Sp  to  denote  the  perturbations 
in  density,  incompressibility  and  rigidity,  respectively.  The  first  Born  approximation 
leads  to 


5u  (x,  t) 


I 


a  *  s  dv 


(25) 


with  (Doornbos  and  Mondt,  1979): 

S  =  -  6pu®  +  (6k--|-6p)  V  (v.u“)  -  VX  7XU® 

+  (v.u®)  v(5K--|-fip)  +  2e® 


(26) 


and  the  strain  tensor  has  the  elements 

•  a,..?) 


Substituting  equations  (2)  and  (3)  for  the  incident  wave  and  the  Green’s  function,  we 
can  write  the  scattered  field  in  a  form  similar  to  equation  (4): 


6u(x,  t)  »  -  /  V  (x)  p  (  a^g^)  dV  (27) 

*  j  =1 

where 
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?!  =  6p/p  ,  =  (v®.v^) 

gj  =  ^k/k  ,  aj  =  -(l--je)  (y^.v”)  (y^.v^) 

gj  =  6p/p  ,  aj  =  |(Y®-v^)  (v®.y^)  +  (y"-Y^)  (v®.v*-) 

-■j  (y®*v“)  (y^.v^)  j  (28) 

Here  y  and  v  are  unit  vectors  in  the  direction  of  wave  propagation  and 
displacement,  respectively,  c  is  the  wave  velocity  (a  for  P,  p  for  S),  and  e=pyc^. 
Quantities  associated  with  the  incident  and  scattered  wave  have  superscripts  °  and 
respectively. 

I^xainining  equations  (27)  and  (4)  reveals  that  the  acoustic  factor 


V.,  =  2c  ds  (v®.v^)  ---  -2  —  (v®.v^) 

c 

is  replaced,  for  elastic  scattering,  by 


(29) 


3 

=  E  (30) 

i  -1 

Thus  we  can  use  equations  (28)-(30)  to  assess  the  validity  of  acoustic  scattering.  T»ie 
velocity  perturbations  can  be  expressed,  assuming  a^fp^=3: 

9  -  Sp  _  5  fiK  4  ftp 

Cl  p  9  1C  Q 

(31) 

-  2-^  =  le.  -  in 

P  P  P 

For  P-*P  scattering,  from  equation  (28)  and  using  cos  0=(y  o  y  >): 

y^j(P)  =  cos<p  -  -i  -  J.  Ocos^ip-l)  At  (32) 

p  9  K  9  \i 


and  for  S-S  scattering: 
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Kl  =  COS(p 


(2cos^<p-l) 


(33) 


Ae.  - 

p 


Ail 

p 


Thus  for  forward  scattering  in  the  so-called  specular  direction,  (p=0  and 
the  Born  approximation,  for  both  P  and  5  waves  (Aki  and  Richards,  1980).  For 
the  error  )  t^cpends  on  the  relative  perturbations  of  p,  k  and  p.  To 

get  a  rough  estimate  of  this  error  we  adopt,  from  Anderson  (1989): 


AA=2A“=i5A£. 

P  a  ■  p 

hence  AlL  =  4.5  —  =1.35  Af. 

p  K  p 


With  these  quantities,  the  relative  error 
e  =  /  ^ao 


increases,  for  P  from  0  for  to  20%  for  0  =  27”,  and  to  100%  for  0=60”. 

The  relative  error  for  S  increases  from  0  for  <p=(P,  to  20%  for  <p  =  18°,  and  to 
100%  for  <p=4(fi.  To  judge  the  significance  of  these  numbers,  consider  a  typical 
tomographic  experiment  with  teleseismic  data  aimed  to  map  upper  mantle  structure 
at  about  100  km  depth  from  data  near  60”epicentral  distance.  Including  half  the 
Fresnel  zone  in  the  diffraction  term  gives  ^  ^  second  period,  and 

*^inax  20  seconds  period.  If  20%  is  taken  as  an  acceptable  relative  error,  then 

the  acoustic  approximation  is  acceptable  for  P  up  to  a  period  of  about  5  seconds. 
The  corresponding  numbers  for  5  are:  ^ma^-lO”  at  1  second  period,  at 

20  seconds  period,  and  the  acceptable  period  range  for  the  acoustic  approximation  is 
up  to  about  4  seconds.  For  longer  periods,  diffraction  tomography  might  proceed 
using  a  prescribed  scaling  relation  between  the  relative  perturbations  of  p,  k  and 
ti¬ 
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A  boundary  perturbation 

Perturbations  of  internal  boundaries  and  of  the  Earth’s  surface  must  be  treated 
different  from  relatively  smooth  velocity  perturbations.  The  difference  arises  both  in 
seismic  tomography  and  in  scattering  theory.  In  tomography  the  travel  time 
perturbation  of  a  wave  being  transmitted  through  or  reflected  from  a  boundary  that  is 
displaced  over  a  distance  h,  is  obtained  by  a  simple  geometrical  analysis  of  the  wave 
front: 


6t  =  iPz-pl)  h 


(34) 


Without  loss  of  generality  we  assume  here  and  in  the  following  that  the  reference 
level  is  horizontal.  Then  li  is  the  vertical  displacement  of  the  boundary  from  the 
reference  level,  and  is  the  vertical  component  of  the  slowness  vector. 

Superscripts  °  and  '  denote  the  incident  and  scattered  waves,  respectively;  in  what 
follows  we  will  use  the  same  notation  to  identify  the  velocities  and  densities  at  the  two 
sides  of  the  boundary.  Using  scattering  theory,  we  want  to  generalize  the  above  result 
in  analogy  to  the  result  (23)  for  velocity  perturbations,  such  that  equation  (34)  still 
gives  the  ray  contribution  in  addition  to  the  diffraction  term. 

In  scattering  theory,  the  generalization  of  reflection/transmission  coefficients  for  a 
plane  boundary  to  coefficients  for  a  rough  boundary  has  been  given  in  the  form  of  a 
perturbation  series  (Doornbos,  1988).  For  our  present  purpose  we  need  the  first  two 
terms  of  this  series,  including  the  Born  approximation,  l^t  the  reference  level  be  the 
z-0  plane,  then  the  boundary  perturbation  is  characterized  by  z=h(^^  where 

are  coordinates  in  the  z=0  plane,  and  the  (non-unit)  normal  to  the  boundary 
is  n  -(-d  hjdx,  -dhldy,  /)T  It  is  appropriate  to  determine  first  a  displacement- 
traction  vector  at  the  boundary.  Let  be  the  displacement  components,  r  ^  the 
associated  stress  tensor,  and  the  modified  traction  components.  Tne 

appropriate  displacement-traction  vector  is,  for  a  solid-solid  interface: 


d  =  {u^,  Uy,  u^.oj  i(ji .  Oy/ i(ji  ,a^/ 


For  a  free  surface: 


d  =  {Ux.UyrU^)’’’ 


where  the  superscript  -  denotes  the  field  below  the  surface.  For  a  solid-liquid 
interface: 
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d=  (u*,Uy,n.u,n.a/i(ti)'^ 


where  are  horizontal  displacement  components  in  the  solid  which  is  here 

assumed  above  the  interface. 

The  solution  for  d  is  obtained  in  the  horizontal  wavenumber  domain.  We  define 


+  00 

D{k^)  =  ff 

—  oo 


(35) 


where  define  an  incident  wave  vector 

A^(k^  )  exp  containing  the  wavenumber  components  of  up-  and  downgoing 

P,  SV  and  SH.  The  complete  solution  for  d(i^,h)  would  require  D(k^  )  for  all 
wavenumbers  k^.  Within  the  present  context  this  is  not  a  practical  solution.  Instead 
we  will  approximate  the  solution  for  d  in  such  a  way  that  it  is  consistent  with  D(k^^) 
for  the  particular  wavenumber  4,'  given  by  the  Green’s  function  (equation  3)  that 
accounts  for  propagation  between  the  scatterer  in  (  and  the  receiver  in  x.  The 
required  approximation  is 


i  ikl)  exp  iiKh)  • 

exp  (36) 


Here  P®,  A®  and  y®  are  reflectivity  matrices  introduced  by  Doornbos  (1988).  They 
play  the  same  role  as  the  matrices  arising  in  the  calculation  of  ordinary  reflection/ 
transmission  coefficients.  In  fact,  P®  is  just  the  matrix  needed  for  calculating  the 
coefficients  for  a  plane  horizontal  boundary.  The  additional  matrices  are  needed  to 
satisfy  conditions  at  the  sloping  boundary.  The  diagonal  matrix  exp(iKh)  contains  the 
vertical  wavenumbers  for  up-  and  downgoing  P,  SV  and  SH.  It  acts  as  a  propagator 
matrix  between  the  reference  level  at  z=0  and  the  boundary  at  z=h. 

Substituting  d(i,  ,h)  in  equation  (35),  expanding  exp(iKh)  in  a  Taylor  series,  and 
retaining  the  zeroth  and  first  order  terms: 


D(k^)  “  D^^Uk^)  +!)«">(*,) 
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it  is  easily  verified  that  the  resulting  expressions  for  and  are  indeed 
equivalent  to  those  given  by  Doornbos  (1988). 

The  scattered  wave  coefficient  vector  B(k^)  can  be  obtained  by  the  matrix  integral 
equation  (c.f.  Doornbos,  1988): 


+  00 

B  (*S)  =  //  exp  {iKh)  I^P^(kl) 

-OO 


-1^  xMJtt) 
ax 


dy 


(37) 


where  P'.  and  y’  are  reflectivity  matrices  in  analogy  to  the  ones  introduced  in 
equation  (36).  We  assume  now  that  there  is  just  one  type  of  incident  wave  A-J>,  such 
that  in  equation  (36)  we  do  not  need  the  complete  inverse  matrix  (P®)  ’(^,®)>  but  only 
the  appropriate  column  (/*),  „'  (k°).  Likewise  we  consider  one  type  of  scattered 
wave  B^,  and  we  need  in  equation  (37)  only  the  appropriate  row  vectors 
and  and  the  appropriate  vertical  wavenumber  10.  We  note  that  for  any  of 

the  up- and  downgoing  waves,  the  vertical  wavenumber  Ar=±<ap^  where  is  the 
vertical  component  of  the  slowness  vector.  In  the  following  we  use  K=o)Z  and 
IO  =  -o)Pj^  (this  follows  from  the  sign  conventions  used  in  Doornbos,  1988). 

Combining  equations  (36)  and  (37)  with  the  above  specifications,  and  expanding 
exp{i(i>Zh)  and  cxp(-ia>p^'li)  in  Taylor  series,  we  get  the  following  result  up  to  order 
one: 


“  f{{Rikl.kl)  >  i(jhRUkl,kt)  +  Vb. J?” (Jtj, *?) ) 

-  OO 

ikl)  exp[-i{kl-kl)-^,] 


(p«) -‘(xn 


l*"/  ^  nO\  -1  / 


RUkl.kl)  =  {{P'‘)-Ukl)  ZP°(kl)  -p]j)(P®)ji  ikt) 


(38) 

(39) 

(40) 


P“(Jt^Jk“) 


P^^ikl)  (P«)-MJtJ) 


x*>(kl)  ' 
,  rVJkJ)  ^ 


i  J 


(P")i'n(*t) 


(41) 
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The  expression  for  the  scattered  field  u(x,t)  requires  a  slight  modification  of  the 
result  (38).  Firstly,  the  coefficient  is  associated  with  the  plane  wave 

component  of  a  Weyl  integral.  It  can  be  associated  with  the  Green’s  tensor  (3)  by 
rewriting,  in  the  frequency  domain: 

0(x,5,w)  =  V  (x)  v^'(0  A^d)  exp(iwrM 

=  -2io>p^c^’pi  V  (x)  Q^(kl)  (i)  exp(ia)r^)  (42) 


where  is  the  appropriate  factor  in  the  Weyl  integral.  Secondly,  the 

wavenumbers  of  the  incident  and  scattered  waves,  and  A,',  vary  with  scatterer 
position  on  the  boundary.  Combining  equations  (38)  and  (42)  and  transforming 
the  result  to  the  time  domain,  we  get 

u,(x,  t)  =  i  W  f  {t-x)  -  hR^f(t~x)]  yl  dS  (43) 

S 


where  W  is  given  by  equation  (6),  and  is  the  vertical  component  of  the  unit 
wave  direction  vector,  i.e.  y//c*.  The  term  with  R  produces  the  zeroth  order 
field.  The  perturbation  is  due  to  the  terms  with  and  /?”,  There  may  be  an 
additional  perturbation  due  to  rapid  lateral  changes  of  R;  we  consider  this  possibility 
in  the  next  section,  but  ignore  the  perturbation  here.  Thus: 


bu^{x,t)  = 


The  term  with  h 
the  integral: 


-/  f/ t-x)  -  vh.R^^  f  (t-x) 
can  be  treated  in  the  same  way  as  equation  (5).  Thus  we  transform 


\ 

y\  dS  (44) 

t 


f  W  h  R^  £{  t-x)  y\  dS  =  -  f  /(t-x)  /  W  R^h  J  dl^  dx 
s  J  L 


(45) 


where  is  the  trajectory  t  =constant  (i.e.  an  isochron,  see  Fig.  2).  Repeating  the 
steps  leading  to  equation  (8),  we  get  here 
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(46) 


-  f  (  t-T 


j  W  R^h  J,  Yi  *  f(t-xj  /  w  R^h  jy,  di,„ 


~  f  W  R^  f  (t-x)  jl  dS 

C  OT 


Note  that  in  the  first  integral,  /?>  is  to  be  evaluated  for  A,*=*o.  From  equations 
(39)  and  (40): 


R^ikl.kt)  =  (pI-pI)  R(kt.kl) 


(47) 


To  evaluate  the  first  integral  of  expression  (46),  we  use  the  result  (10)  for  an 
infinitesimal  element  of  surface  normal  to  the  scattered  ray.  The  factor  is 
needed  for  the  transformation  to  an  element  of  surface  on  the  reference  plane.  Thus 
(c.f.  equation  11): 


dg!  Yz 


2il 

■icTl/2 


(48) 


Noting  that  the  spreading  matrix  g®  in  equation  (12)  applies  to  the  incident  wave, 
we  have  to  multiply  (g”l  by  the  factor  upon  transmission  through  the 

boundary  (Cerveny,  1985).  The  amplitude  factor  A  in  equation  (16)  is  similarly 
modified  by  the  factor 

Thus  using  equation  (48),  (12),  (13),  (14)  and  (16),  and  substituting  (47), 


6u  “(X,  t) 


/  (t  /  w  R^  h  J„y\  dl,„ 
(p°-pi)  -b  u  "(x,  t) 


(49) 


It  is  argued,  as  in  the  previous  section,  that  the  second  term  of  expression  (46) 
represents  diffraction  from  the  btiundary  of  surface  5,  and  this  effect  should  be 
deleted  in  practical  circumstances  (the  implicit  assumption  being  h=constant  outside 
5).  To  evaluate  the  last  integral  of  (46),  we  repeat  the  development  summarized  by 
equations  (17)-(20),  except  that  we  use  2-D  instead  of  3-D  vectors,  i.e.  we  express  our 
results  in  terms  of  horizontal  slowness  vectors  and  />,>,  where  <op°=k^^  and 


16 


and  an  area  bounded  by  St  =coiLsiant  close  to  the  ray  can  be  excluded  from 
integration. 


Summarizing  the  results  expressed  by  (45),  (46)  and  (50),  we  can  rewrite  equation 
(44): 


6u^(x,t)  -  Su/Cx,  t)  +  6uf(x,t) 
=  -  h  d™(x,  t) 


-/ 


W  vh  . 


'  (pI-pI) 

.  |Pe“-pJP 


f  it-x)  Yz  dS 


(51) 


In  the  frequency  domain: 


bU^{x,(j))  -  bUi{x,(o)  +5£//(x, <i)) 


=  io)  (p°-pl)  h  u"(x,o)) 


+  i<i}  F(a))  f  W  vh 
S 


exp  (i(i)T)  Y^  dS  (52) 


The  Rytov  approximation  is: 


In  -  io>  [(p°-pj)  h 

C//(x,w)  \ 


Vj(x)A(x)  s 


/  6^ 


vh 


(P^P^ 

|p2-PeN^ 


exp  [io)  (t-tJ  ]  Y^  dS 


(53) 


Note  the  similarity  between  equations  (51)-(53),  and  equations  (21)-(23)  for  acoustic 
scattering.  The  first  term  in  these  equations  is  the  usual  delay  time  of  seismic 
tomography.  However,  the  diffraction  integral  for  a  boundary  perturbation  contains 
an  additional  term  v/i.lf"  that  has  no  counterpart  in  acoustic  scattering.  This  term 


17 


is  needed  to  satisfy  conditions  at  a  non-horizontal  boundary.  Examination  of  and 
(equations  40  and  41)  suggests  that  the  two  diffraction  terms  can  be  of 
comparable  magnitude. 

All  of  the  remarks  regarding  the  implementation  of  acoustic  diffraction  tomography 
following  equation  (23)  are  also  relevant  to  boundary  tomography  following  equation 
(53).  Thus  the  integration  surface  5  can  be  chosen  to  correspond  to  (an  appropriate 
fraction  of)  the  Fresnel  zone,  and  the  time  r  may  have  to  be  calculated  iteratively 
using  previous  tomographic  results  for  h.  The  calculation  of  /?•  and  7?"  in 
equation  (53)  is  more  elaborate  than  the  calculation  of  the  corresponding  factor 
(yO.  v’)  for  acoustic  scattering.  On  the  other  hand,  and  are  frequency 
independent,  and  the  integration  in  equation  (53)  is  2-D  rather  than  the  3-D 
diffraction  integral  in  equation  (23). 
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Perturbation  of  boundary  conditions 


We  reconsider  the  fi'^st  term  of  the  scattered  field,  in  equation  (43).  The  integral  with 
the  factor  R  can  he  evaluated  in  the  same  way  as  the  integral  with  R\  using 
equations  (45)-(50).  This  lettds  to: 


f  W  R  f  (t-x)  Yz  =  u/’(x,  t)  hui(x.  t) 
S 


and 


(Pc  pi ) 

Ipi-pil^ 


V  R  f(t-x)  dS 


(54) 


(55) 


In  the  frequency  domain: 


butix.ui)  -  F(a))  /  w  ■ 

5  \pI-pI\'^ 


V  R  oxp  (  /(i)t )  dS 


(56) 


In  order  to  assess  the  significance  of  due  to  vR,  relative  to  due  to  v/i 

in  equation  (52),  we  write 


vR  =  ftp*  +  ftp 


8p*  8p 


5k‘  +  6k 


8k* 


8k 


+  ftp* 


8p*  8p 


♦  ^  5m- 


(57) 


Hence  perturbations  Sp,  Sk,  Sp  along  the  boundary  produce  variations  6R  of  the 
order 


kAZ  ,  R.P±  ,  rM 

P  K  p 


On  the  other  hand,  a  comparable  variation  due  to  a  perturbation  6h  can  be 
deduced  from  equation  (52),  and  is  of  the  order 


b)R  (p°-pi)  6/3 
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Consider  for  example  a  bottomside  reflection  (PKKP)  at  the  core-mantle  boundary. 
A  typical  ray  parameter  value  might  be  2.5  s/d.  Take  Then 


6h  “  1.5  !)h  at  IHz. 


Thus  relative  perturbations  6p/p,  6kJk,  S^/p  of  the  order  of  10%  would  be  needed 
to  simulate  the  effect  of  a  boundary  perturbation  Sit  of  the  order  of  100  m.  This 
i)rder  of  magnitude  argument  suggests  that  in  the  above  situation,  the  effect  of  a 
boundary  perturation  Sli  dominates  that  of  moderate  perturbations  of  the  elastic 
constants  and  density  along  the  boundary.  However  the  diffraction  term  due  to  h 
and  6h  decreases  with  increasing  wave  period  and  incident  angle;  hence  in  practice 
the  relative  importance  of  SR  and  Sli  would  have  to  be  assessed  for  the  actual 
situation  at  hand. 
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Conclusions 


Diffraction  tomography  can  be  formulated  in  such  a  way  that  it  makes  explicit  the 
phase  delay  perturbation  from  ray  theory  which  is  the  basis  of  ordinary  seismic 
tomography.  The  phase  delay  term  depends  on  the  velocity  perturbations  along  the 
ray  in  both  acoustc  and  elastic  media,  and  when  the  ray  crosses  a  boundary  an 
additional  phase  delay  is  induced  that  is  proportional  to  the  boundary  level 
perturbation. 

The  additional  diffraction  term  involves  both  phase  and  amplitude  perturbations.  The 
diffraction  depends  on  the  gradients  of  the  velocity  perturbation  in  an  acoustic 
medium,  the  gradients  of  the  elastic  and  density  perturbations  in  an  elastic  medium, 
and  the  gradients  of  the  boundary  perturbations  the  wave  is  crossing.  The  diffraction 
term  arises  also  in  circumstances  when  the  primary  wave  is  cut  off  from  the  receiver, 
for  example  in  a  shadow  zone.  The  formulas  for  the  diffraction  term  due  to  acoustic, 
elastic  and  boundary  perturbations  are  very  similar,  but  the  boundary  perturbation 
'■equires  an  additional  term  to  satisfy  conditions  at  a  sloping  boundary.  In  the  present 
formulation,  numerical  diffraction  from  the  nonphysical  boundary  of  the  region  under 
study  appears  as  a  separate  term  and  can  thus  be  easily  removed.  This  is  consistent 
w'ith  the  assumption  of  constant  perturbations  outside  this  region,  in  contrast  to  the 
conventional  formulation  which  implies  that  the  perturbation  drops  to  zero  outside 
the  region  under  study.  This  aspect  is  especially  important  when  inverting  travel 
times,  or  relatively  short  waveform  sections. 

The  difference  between  an  elastic  and  an  acoustic  medium  is  immaterial  in  ordinary 
seismic  tomography  since  only  velocity  perturbations  can  be  retrieved.  In  diffraction 
tomography,  the  difference  may  be  neglected  when  short-period  data  are  used.  In  a 
realistic  geometry  that  is  representative  of  tomography  for  upper  mantle  strucutre,  the 
acceptable  period  range  was  found  to  be  about  5  seconds  for  P  and  about  4  seconds 
for  5,  if  half  the  Fresnel  zone  is  included  and  the  acceptable  relative  error  is  20  %. 
For  longer  periods  a  preferred  alternative  would  be  to  invert  for  scaled  perturbations 
of  the  elastic  constants  and  density. 

Diffraction  from  a  boundary  can  be  induced  not  only  by  gradients  of  the  boundray 
level,  but  also  by  lateral  gradients  of  the  elastic  constants  and  densities  along  the 
boundary.  The  relative  importance  of  such  perturbations  of  boundary  conditions 
depends  on  the  wave  period  and  incidence  angle.  A  rough  order  of  magnitude 
calculation  suggests  that  in  some  recent  tomographic  studies  of  the  core-mantle 
boundary,  the  effect  of  boundary  level  perturbations  is  probably  more  important  than 
the  effect  of  up  to  moderate  perturbations  of  the  elastic  constants  and  densities. 
However  one  might  have  to  rea.ssess  this  conclusion  in  other  circumstances. 
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Fig.  1  -  Schematic  diagram  showing  the  ray  between  0  (source)  and  1  (receiver),  with 
stationary  travel  time  t^.  Also  shown  is  a  section  of  isochron  surface  S^,  and 
the  surface  normal  to  the  ray  is  S^.  Surface  S^  is  defined  by  all  points  ( 
such  that  the  total  travel  time  of  the  two  rays  connecting  (f,  0)  and  {(,  1)  is 

T. 
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Fig.  2  -  Schematic  diagram  showing  the  reflected  ray  between  0  (source)  and  1 

(receiver),  with  stationary  travel  time  Also  shown  is  the  isochron  on 
the  reflecting  surface  5,  and  the  surface  normal  to  the  reflected  ray  is  5^.  A 
similar  diagram  applies  to  a  refracted  ray. 
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